﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "euler2")]
    public static unsafe void euler2(double t, double h, IntPtr y_ptr, IntPtr z_ptr, double eps, IntPtr f_abc_ptr)
    {
        double* y = (double*)y_ptr.ToPointer();
        double* z = (double*)z_ptr.ToPointer();
        f_abc = Marshal.GetDelegateForFunctionPointer<delegatefunc_abc>(f_abc_ptr);

        euler2(t, h, y, z, eps);
    }

    /// <summary>
    /// 求解二阶初值连分式法
    /// 改进欧拉公式以h为步长积分m步
    /// f计算二阶微分方程的右端函数f(t,y,z)。
    /// </summary>
    /// <param name="t">自变量起点值</param>
    /// <param name="h">步长</param>
    /// <param name="y">存放函数初值。返回终点函数值。</param>
    /// <param name="z">存放函数一阶导数初值。返回终点函数一阶导数值。</param>
    /// <param name="m">步数</param>
    public static unsafe void euler21(double t, double h, double* y, double* z, int m)
    {
        int j;
        double x, yy, zz, yc, zc, yk1, yk2, zk1, zk2;
        yy = *y; zz = *z;
        for (j = 0; j <= m - 1; j++)
        {
            x = t + j * h;
            //计算yK1
            yk1 = zz;
            //计算zK1
            zk1 = f_abc(x, yy, zz);
            x = t + (j + 1) * h;
            //预报t[j+1]处的y值
            yc = yy + h * zk1;
            //预报t[j+1]处的z值
            zc = zz + h * yk1;
            //计算yK2
            yk2 = zc;
            //计算zk2
            zk2 = f_abc(x, yc, zc);
            //计算t[j+1]处的y值
            yy = yy + h * (yk1 + yk2) / 2.0;
            //计算t[j+1]处的z值
            zz = zz + h * (zk1 + zk2) / 2.0;
        }
        *y = yy; *z = zz;
        return;
    }

#if __DUPLICATED__
    // 求解二阶初值Euler方法.cpp
    // 改进欧拉公式以h为步长积分m步
    // t自变量起点值
    // h步长
    // y存放函数初值。返回终点函数值。
    // z存放函数一阶导数初值。返回终点函数一阶导数值。
    // m步数
    // f二阶微分方程右端函数f(t,y,z)。
    //  void euler21(double t,double h,double *y,double *z,int m,	         double (*f)(double,double,double))
    public static unsafe void euler21(double t, double h, double* y, double* z, int m)
    {
        int j;
        double x, yy, zz, yc, zc, yk1, yk2, zk1, zk2;

        yy = *y; zz = *z;
        for (j = 0; j <= m - 1; j++)
        {
            x = t + j * h;
            //计算yK1
            yk1 = zz;
            //计算zK1
            zk1 = (*f)(x, yy, zz);
            x = t + (j + 1) * h;
            //预报t[j+1]处的y值
            yc = yy + h * zk1;
            //预报t[j+1]处的z值
            zc = zz + h * yk1;
            //计算yK2
            yk2 = zc;
            //计算zk2
            zk2 = (*f)(x, yc, zc);
            //计算t[j+1]处的y值
            yy = yy + h * (yk1 + yk2) / 2.0;
            //计算t[j+1]处的z值
            zz = zz + h * (zk1 + zk2) / 2.0;
        }
        *y = yy; *z = zz;
        return;
    }
#endif

    // Euler方法变步长积分一步二阶初值问题
    // t自变量起点值
    // h步长
    // y存放函数初值。返回终点函数值。
    // z存放函数一阶导数初值。返回终点函数一阶导数值。
    // eps精度要求
    // f二阶微分方程右端函数f(t,y,z)。
    public static unsafe void euler2(double t, double h, double* y, double* z, double eps)
    {
        int m;
        double p, ya, za, yb, zb;
        m = 1;
        p = 1.0 + eps;
        ya = *y; za = *z;
        //跨一步计算
        euler21(t, h, &ya, &za, m);
        while (p > eps)
        {
            yb = *y; zb = *z; m = m + m; h = h / 2.0;
            //跨m步计算
            euler21(t, h, &yb, &zb, m);
            //取误差
            p = Math.Abs(yb - ya);
            za = zb; ya = yb;
        }
        *y = ya; *z = za;
        return;
    }

    /*
    // 求解二阶初值Euler方法例1
      int main()
      { 
          int j;
          double t,h,eps,y,z;
          double  euler2_f(double, double, double);
          y=0.0; z=0.701836;
          t=0.0; h=0.1; eps=0.0000001;
          cout <<"t = " <<setw(6) <<t;
          cout <<setw(6) <<"y = " <<setw(10) <<y;
          cout <<setw(6) <<"z = " <<setw(10) <<z;
          cout <<endl;
          for (j=1; j<=10; j++)
          { 
              euler2(t,h,&y,&z,eps,euler2_f);
              t=t+h;
              cout <<"t = " <<setw(6) <<t;
              cout <<setw(6) <<"y = " <<setw(10) <<y;
              cout <<setw(6) <<"z = " <<setw(10) <<z;
              cout <<endl;
          }
          return 0;
      }
    // 计算二阶微分方程右端函数f(t,y,z)
      double euler2_f(double t, double y, double z)
      { 
          double d;
          d=t+y;
          return(d);
      }
    */
    /*
    // 求解二阶初值Euler方法例2
      int main()
      { 
          int j;
          double t,h,eps,y,z;
          double  euler2_f(double, double, double);
          y=1.0; z=0.0;
          t=0.0; h=0.1; eps=0.0000001;
          cout <<"t = " <<setw(6) <<t;
          cout <<setw(6) <<"y = " <<setw(10) <<y;
          cout <<setw(6) <<"z = " <<setw(10) <<z;
          cout <<endl;
          for (j=1; j<=10; j++)
          { 
              euler2(t,h,&y,&z,eps,euler2_f);
              t=t+h;
              cout <<"t = " <<setw(6) <<t;
              cout <<setw(6) <<"y = " <<setw(10) <<y;
              cout <<setw(6) <<"z = " <<setw(10) <<z;
              cout <<endl;
          }
          return 0;
      }
    // 计算二阶微分方程右端函数f(t,y,z)
      double euler2_f(double t, double y, double z)
      { 
          double d;
          d = (6*t-3.0+t*z+3*y)/(1.0+t*t);
          return(d);
      }
    */
}

